Figures and supplemental figures for ecological resilience ms

Author

Botsford, Hastings, White, and Kilduff

List of tables in “data/sim_out.duckdb”:

Config and set up

This .rmd file can be run as a script that produces the main text figues and the supporting figures for the manuscript Marine reserves can provide buffering against environmental fluctuations to overexploited but not sustainably harvested fisheries.

In RStudio, click on Code -> Run Region –> Run All to run the script and create all the files.

To render an .html file from RStudio, click the Render button.

To export to pdf, follow the instructions here to install a Quarto “friendly” Tex distribution. Then click the small down arrow of the Render button and select Render as pdf

Exporting to a Word (.docx) file requires that MS Word/Office is installed, then click the small down arrow of the Render button and select Render as MS Word .

Main text figures

# Set table names to be used to make the main figures for blue rockfish 


table_name_brf_nd_mei <- "no_dispersal_blue_rockfish_mei_135_75"
table_name_brf_lp_mei <- "larval_pool_blue_rockfish_mei_135_75"
table_name_brf_nd_white <- "no_dispersal_blue_rockfish_white_135_75"
table_name_brf_lp_white <- "larval_pool_blue_rockfish_white_135_75"

# Get meta data descriptins from the table names and the pretty 
# versions for plotting

brf_nd_mei_current <- get_current_exp(table_name = table_name_brf_nd_mei)
brf_nd_mei_nice <- get_nice_exp(table_name = table_name_brf_nd_mei)

brf_nd_white_current <- get_current_exp(table_name = table_name_brf_nd_white)
brf_nd_white_nice <- get_nice_exp(table_name = table_name_brf_nd_white)

brf_lp_mei_current <- get_current_exp(table_name = table_name_brf_lp_mei)
brf_lp_mei_nice <- get_nice_exp(table_name = table_name_brf_lp_mei)

brf_lp_white_current <- get_current_exp(table_name = table_name_brf_lp_white)
brf_lp_white_nice <- get_nice_exp(table_name = table_name_brf_lp_white)

## Blue rockfish reproductive biology and fishery params

# use species specific vB params to create a weight at age data.frame that can be used to convert abundance at age to biomass at age
species_name <- brf_nd_mei_nice[[2]]
brf_waa_g <- sim_species_parms[[species_name]][["biom_const"]] * sim_species_derived_vars[[species_name]][["length_at_age"]] ^
                    sim_species_parms[[species_name]][["biom_exp"]]

brf_waa_g_df <- data.frame(age = 1: sim_species_parms[[species_name]]$A_max,
                           waa_g = brf_waa_g)

# use functions to create useful data structures
# FLEP to F reference tables by species
spec_fleps_2_fs_brf <-
  get_fs_from_fleps(species = species_name,
                    flep_inds = species_flep_fished_ind,
                    f_vals =  fishing_mortality_values)

Figure 3

Contour plots showing the buffering effects of harvest rate and marine reserves on the abundance of a fished population, when there is no larval dispersal between the fished and reserve habitat patches.

Species: Blue rockfish

Dispersal: No Dispersal

Environmental noise: ENSO (MEI)

sim_res_out_brf_nd_mei <-
  qry_get_sim_data_for_analysis(table_name = table_name_brf_nd_mei)
# sim_res_out_brf_nd_mei

F_brf_maxes_yield_brf_nd_mei <- get_f_maxes_yield(input_df = sim_res_out_brf_nd_mei,
                                         fleps_2_fs = spec_fleps_2_fs_brf)
# F_brf_maxes_yield_brf_nd_mei

out_fs_brf_nd_mei <- df_add_fs_out(
  input_df = sim_res_out_brf_nd_mei,
  spec_fleps_2_fs = spec_fleps_2_fs_brf,
  F_maxes_yield = F_brf_maxes_yield_brf_nd_mei,
  waa_g_df = brf_waa_g_df
)
# out_fs_brf_nd_mei

med_yield_brf <- get_median_yield(input_df = out_fs_brf_nd_mei) 
# med_yield_brf

med_biomass_brf <- get_median_biomass(input_df = out_fs_brf_nd_mei)
# med_biomass__brf_nd_mei

out_yld_bm_sims_f_frac_yr_brf_nd_mei <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = out_fs_brf_nd_mei)
# out_yld_bm_sims_f_frac_yr_brf_nd_mei

sum_stats_f_frac_df_brf_nd_mei <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = out_yld_bm_sims_f_frac_yr_brf_nd_mei,
    median_yield = med_yield_brf,
    median_biomass = med_biomass_brf
  )
# sum_stats_f_frac_df_brf_nd_mei

## Plotting

p1 <-
  plot_contour_bm_below_thresh(sum_stats_f_frac_df = sum_stats_f_frac_df_brf_nd_mei,
                               spec_disp_noise = brf_nd_mei_nice)

p2 <-
  plot_contour_yld_below_thresh(sum_stats_f_frac_df = sum_stats_f_frac_df_brf_nd_mei,
                                spec_disp_noise = brf_nd_mei_nice)

fig3 <-
  plot_grid(
    p1,
    p2,
    labels = letters[1:2],
    label_size = 9,
    byrow = TRUE,
    nrow = 2,
    ncol = 1
  )

print(fig3)

ggsave(filename = "output/main_figs/fig3-brf-nd-mei-biomass-yield.png",
       plot = fig3,
       height = 10,
       width = 8)
ggsave(filename = "output/main_figs/fig3-brf-nd-mei-biomass-yield.svg",
       plot = fig3,
       height = 10,
       width = 8)

Figure 4

Contour plots showing the buffering effects of harvest rate and marine reserves on the abundance of a fished population, when there is larval pool dispersal between the fished and reserve habitat patches.

Species: Blue rockfish

Dispersal: Larval Pool Dispersal

Environmental noise: ENSO (MEI)

## Qry 

sim_res_out_brf_lp_mei <-
  qry_get_sim_data_for_analysis(table_name = table_name_brf_lp_mei)

F_brf_maxes_yield_brf_lp_mei <- get_f_maxes_yield(input_df = sim_res_out_brf_lp_mei,
                                         fleps_2_fs = spec_fleps_2_fs_brf)
# F_brf_maxes_yield_brf_lp_mei

out_fs_brf_lp_mei <- df_add_fs_out (input_df = sim_res_out_brf_lp_mei,
                         spec_fleps_2_fs = spec_fleps_2_fs_brf,
                         F_maxes_yield = F_brf_maxes_yield_brf_lp_mei, 
                         waa_g_df = brf_waa_g_df)
# out_fs_brf_lp_mei

out_yld_bm_sims_f_frac_yr_brf_lp_mei <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = out_fs_brf_lp_mei)
# out_yld_bm_sims_f_frac_yr_brf_lp_mei

sum_stats_f_frac_df_brf_lp_mei <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = out_yld_bm_sims_f_frac_yr_brf_lp_mei,
    median_yield = med_yield_brf,
    median_biomass = med_biomass_brf
  )
# sum_stats_f_frac_df_brf_lp_mei

## Plotting

p1 <-
  plot_contour_bm_below_thresh(sum_stats_f_frac_df = sum_stats_f_frac_df_brf_lp_mei,
                               spec_disp_noise = brf_lp_mei_nice)

p2 <-
  plot_contour_yld_below_thresh(sum_stats_f_frac_df = sum_stats_f_frac_df_brf_lp_mei,
                                spec_disp_noise = brf_lp_mei_nice)

fig4 <-
  plot_grid(
    p1,
    p2,
    labels = letters[1:2],
    label_size = 9,
    byrow = TRUE,
    nrow = 2,
    ncol = 1
  )

print(fig4)

ggsave(
  filename = "output/main_figs/fig4-brf-lp-mei-biomass-yield.png",
  plot = fig4,
  height = 10,
  width = 8
)
ggsave(
  filename = "output/main_figs/fig4-brf-lp-mei-biomass-yield.svg",
  plot = fig4,
  height = 10,
  width = 8
)

Figure 5 and 6

Representative timeseries of total population biomass and yield ( under different combinations of fishing rate and marine reserve protection.

Species: Blue rockfish

Dispersal: Larval pool

Environmental noise: ENSO (MEI)

biom_ts_plt <- facet_plot_bm_ts_plots(
  bio_yield_df = out_fs_brf_lp_mei,
  biomass_thresh = med_biomass_brf,
  yield_thresh = med_yield_brf,
  exp_deets = brf_lp_mei_nice,
  r_seed = 34,
  frac_reserves = c(0, 0.15, 0.30),
  f_fmsys = c(0, 0.5, 1.00, 2.)
)

yield_ts_plt <- facet_plot_yld_ts_plots(
  bio_yield_df = out_fs_brf_lp_mei,
  biomass_thresh = med_biomass_brf,
  yield_thresh = med_yield_brf,
  exp_deets = brf_lp_mei_nice,
  r_seed = 34,
  frac_reserves = c(0, 0.15, 0.30),
  f_fmsys = c(0.5, 1.00, 2.)
)

print(yield_ts_plt)

print(biom_ts_plt)

ggsave(filename = "output/main_figs/fig5-biomass-timeseries-brf-lp-mei.png", plot = biom_ts_plt,
       height = 8, width = 10)
ggsave(filename = "output/main_figs/fig5-biomass-timeseries-brf-lp-mei.svg", plot = biom_ts_plt,
       height = 8, width = 10)

ggsave(filename = "output/main_figs/fig6-yield-timeseries-brf-lp-mei.png", plot = yield_ts_plt,
       height = 8, width = 10)
ggsave(filename = "output/main_figs/fig6-yield-timeseries-brf-lp-mei.svg", plot = yield_ts_plt,
       height = 8, width = 10)

Figure 5 with text

### Biomass

frac_reserves <- c(0, 0.15, 0.3)
f_fmsys <- c(0, 0.5, 1.00, 2)

u_f_fmsys <- unique(sum_stats_f_frac_df_brf_lp_mei[[2]]$f_fmsy)
  
  u_f_fmsys <-
    purrr::map_dbl(f_fmsys,  ~ u_f_fmsys[which.min(abs(u_f_fmsys - .x))])

summary_stats_df_plot_bm <- sum_stats_f_frac_df_brf_lp_mei[[2]] %>%
  dplyr::filter(reduce(map(
    frac_reserves, near, x = as.numeric(reserve_frac), tol = 1e-4
  ), `|`)) %>%
  dplyr::filter(reduce(map(
    u_f_fmsys, near, x = f_fmsy, tol = 1e-4
  ), `|`))

summary_stats_df_plot_yld <- sum_stats_f_frac_df_brf_lp_mei[[2]] %>%
  dplyr::filter(reduce(map(
    frac_reserves, near, x = as.numeric(reserve_frac), tol = 1e-4
  ), `|`)) %>%
  dplyr::filter(reduce(map(
    u_f_fmsys[-1], near, x = f_fmsy, tol = 1e-4
  ), `|`))

# Location to place text
x_text_val <- 135/2

biom_ts_plt <- biom_ts_plt +
  geom_text(
    data = summary_stats_df_plot_bm,
    x = x_text_val,
    y = 3.4,
    size = 3,
    aes(label = mean_rel_biomass_lab)
  ) +
  geom_text(
    data = summary_stats_df_plot_bm,
    x = x_text_val,
    y = 3.1,
    size = 3,
    aes(label = sd_rel_biomass_lab)
  ) +
  geom_text(
    data = summary_stats_df_plot_bm,
    x = x_text_val,
    y = 2.8,
    size = 3,
    aes(label = cv_rel_biomass_lab)
  ) +
  geom_text(
    data = summary_stats_df_plot_bm,
    x = x_text_val,
    y = 2.5,
    size = 3,
    aes(label = yrs_below_biomass_thresh_lab)
  )


print(biom_ts_plt)

ggsave(
  filename = "output/main_figs/fig5-biomass-timeseries-brf-lp-mei-text.png",
  plot = biom_ts_plt,
  height = 8,
  width = 14
)
ggsave(
  filename = "output/main_figs/fig5-biomass-timeseries-brf-lp-mei-text.svg",
  plot = biom_ts_plt,
  height = 8,
  width = 14
)

Figure 6 with text

# Yield

yield_ts_plt <- yield_ts_plt +
  geom_text(
    data = summary_stats_df_plot_yld,
    x = x_text_val,
    y = 1.9,
    size = 2,
    aes(label = mean_rel_yield_lab)
  ) +
  geom_text(
    data = summary_stats_df_plot_yld,
    x = x_text_val,
    y = 1.8,
    size = 2,
    aes(label = sd_rel_yield_lab)
  ) +
  geom_text(
    data = summary_stats_df_plot_yld,
    x = x_text_val,
    y = 1.7,
    size = 2,
    aes(label = cv_rel_yield_lab)
  ) +
  geom_text(
    data = summary_stats_df_plot_yld,
    x = x_text_val,
    y = 1.6,
    size = 2,
    aes(label = yrs_below_yield_thresh_lab)
  )

print(yield_ts_plt)

ggsave(
  filename = "output/main_figs/fig6-yield-timeseries-brf-lp-mei-text.png",
  plot = yield_ts_plt,
  height = 8,
  width = 14
)
ggsave(
  filename = "output/main_figs/fig6-yield-timeseries-brf-lp-mei-text.svg",
  plot = yield_ts_plt,
  height = 8,
  width = 14
)

Figure 7

Representative timeseries of fishery yield (black curves) in the same harvest and marine reserve scenarios depicted in Fig. 5.

Species: Blue rockfish

Dispersal: Larval pool

Environmental noise: ENSO (MEI)

t_c <- species_parms[which(species_parms$species == species_name),]$A_fish

u_f_fmsys <- unique(out_fs_brf_lp_mei$f_fmsy)

f_fmsys_targs <- c(0.5, 1.00, 2.)

u_f_fmsys <-
  purrr::map_dbl(f_fmsys_targs,  ~ u_f_fmsys[which.min(abs(u_f_fmsys - .x))])

recs_spec_disp_noise_ffmsy_plt <- out_fs_brf_lp_mei %>%
  dplyr::filter(sim_num == 1, age == 1) %>%
  dplyr::mutate(reserve_frac = as.numeric(reserve_frac)) %>%
  dplyr::filter(reduce(map(
    frac_reserves, near, x = reserve_frac, tol = 1e-4
  ), `|`),
  reduce(map(
    u_f_fmsys, near, x = f_fmsy, tol = 1e-4
  ), `|`)) %>%
  mutate(year = year - 500,
         lagged_year = year + t_c) 


yield_spec_disp_noise_prep_plt <- out_fs_brf_lp_mei %>% 
  dplyr::filter(sim_num == 1) %>%
  dplyr::mutate(reserve_frac = as.numeric(reserve_frac)) %>%
  dplyr::group_by(reserve_frac, f_fmsy, year) %>%
  dplyr::summarize(annual_yield = sum(yield)) %>%
  dplyr::ungroup() %>%
  dplyr::filter(reduce(map(
    frac_reserves, near, x = reserve_frac, tol = 1e-4
  ), `|`),
  reduce(map(
    u_f_fmsys, near, x = f_fmsy, tol = 1e-4
  ), `|`))  %>%
  mutate(year = year - 500)

plt_start_yr <- 100
plt_end_yr <- 135

scaled_yield_plt <- yield_spec_disp_noise_prep_plt %>%
  group_by(reserve_frac, f_fmsy) %>%
  mutate(scaled_yield = scale(annual_yield)[,1]) %>%
  ungroup()

scaled_recs_plt <- recs_spec_disp_noise_ffmsy_plt %>% 
  group_by(reserve_frac, f_fmsy) %>%
  mutate(scaled_recruits = scale(biomass_at_age)[,1]) %>%
  ungroup()

scaled_yield_plt_4corr <- scaled_yield_plt %>%
  mutate(reserve_frac = as.character(reserve_frac), 
         f_fmsy = as.character(f_fmsy)) %>%
  select(reserve_frac, f_fmsy, year, scaled_yield)

scaled_recs_plt_4corr <- scaled_recs_plt %>%
  mutate(reserve_frac = as.character(reserve_frac), 
         f_fmsy = as.character(f_fmsy))  %>%
  select(reserve_frac, f_fmsy, lagged_year, scaled_recruits)

corr_vals <- scaled_yield_plt_4corr %>%
  left_join(scaled_recs_plt_4corr,
            by = c("reserve_frac" = "reserve_frac",
                   "f_fmsy" = "f_fmsy",
                   "year" = "lagged_year")) %>%
  filter(.,complete.cases(reserve_frac, f_fmsy, year, scaled_yield, scaled_recruits)) %>%
  group_by(reserve_frac, f_fmsy) %>%
  summarise(corr = cor(scaled_yield, scaled_recruits))

lagged_recs_yield_plt <- ggplot(scaled_yield_plt %>% 
  filter(year > plt_start_yr, 
         year <= plt_end_yr), 
       aes(x = year, y = scaled_yield)) + 
  facet_grid(reserve_frac ~ f_fmsy,
               labeller = labeller(f_fmsy = label_both,
                                   reserve_frac = label_both)) +
  geom_line(linewidth = 1.2) +
  geom_line(data = scaled_recs_plt %>% filter(lagged_year > plt_start_yr, 
         lagged_year <= plt_end_yr), aes(x = lagged_year, y = scaled_recruits), 
            colour = "blue", linewidth = 1.2) +
  geom_text(
      data = corr_vals,
      x = 1.1 * plt_start_yr ,
      y = 3,
      aes(label = paste0("Corr: ", round(corr, 3)))
  ) +
  ylab("Scaled yield (black) and lagged + scaled recruits (blue)") +
  xlab("Year") +
  theme_bw()

print(lagged_recs_yield_plt)

ggsave(
  filename = "output/main_figs/fig7-lagged-recs-yield-timeseries-lw-brf-lp-mei.png",
  plot = lagged_recs_yield_plt,
  height = 8,
  width = 10
)

ggsave(
  filename = "output/main_figs/fig7-lagged-recs-yield-timeseries-lw-brf-lp-mei.svg",
  plot = lagged_recs_yield_plt,
  height = 8,
  width = 10
)

Supplemental figure prep

table_name_chr_nd_white <- "no_dispersal_china_rockfish_white_135_75"
table_name_cab_nd_white <- "no_dispersal_cabezon_white_135_75"
table_name_chr_lp_white <- "larval_pool_china_rockfish_white_135_75"
table_name_cab_lp_white <- "larval_pool_cabezon_white_135_75"

table_name_chr_nd_mei <- "no_dispersal_china_rockfish_mei_135_75"
table_name_cab_nd_mei <- "no_dispersal_cabezon_mei_135_75"
table_name_chr_lp_mei <- "larval_pool_china_rockfish_mei_135_75"
table_name_cab_lp_mei <- "larval_pool_cabezon_mei_135_75"


# Get meta data descriptins from the table names and the pretty 
# versions for plotting

chr_nd_white_nice <- get_nice_exp(table_name = table_name_chr_nd_white)
cab_nd_white_nice <- get_nice_exp(table_name = table_name_cab_nd_white)
chr_lp_white_nice <- get_nice_exp(table_name = table_name_chr_lp_white)
cab_lp_white_nice <- get_nice_exp(table_name = table_name_cab_lp_white)
chr_nd_mei_nice <- get_nice_exp(table_name = table_name_chr_nd_mei)
cab_nd_mei_nice <- get_nice_exp(table_name = table_name_cab_nd_mei)
chr_lp_mei_nice <- get_nice_exp(table_name = table_name_chr_lp_mei)
cab_lp_mei_nice <- get_nice_exp(table_name = table_name_cab_lp_mei)


## Qry 

run_qry_get_sims <- function(table_name, burn_in = 500) {
  con <- DBI::dbConnect(duckdb::duckdb(),
                        dbdir = "data/sim_out.duckdb",
                        read_only = FALSE)
  qry <-
    paste0(
      "SELECT  sim_num, reserve_frac, flep_ratio, year, age,
              sum(number) as abundance,
              sum(yield) as yield,
              FROM '",
      table_name,
      "' WHERE year > 500
              GROUP BY sim_num, reserve_frac, flep_ratio, year, age
              ORDER BY sim_num, reserve_frac, flep_ratio, year, age
              "
    )
  
  sim_res_out <- dbGetQuery(conn = con, qry)
  sim_res_out$year = sim_res_out$year - 500
  DBI::dbDisconnect(conn = con, shutdown=TRUE)
  return(sim_res_out)
}


chr_nd_white_df <-
  run_qry_get_sims(table_name = table_name_chr_nd_white, burn_in = 500)
cab_nd_white_df <-
  run_qry_get_sims(table_name = table_name_cab_nd_white, burn_in = 500)
chr_lp_white_df <-
  run_qry_get_sims(table_name = table_name_chr_lp_white, burn_in = 500)
cab_lp_white_df <-
  run_qry_get_sims(table_name = table_name_cab_lp_white, burn_in = 500)
chr_nd_mei_df <-
  run_qry_get_sims(table_name = table_name_chr_nd_mei, burn_in = 500)
cab_nd_mei_df <-
  run_qry_get_sims(table_name = table_name_cab_nd_mei, burn_in = 500)
chr_lp_mei_df <-
  run_qry_get_sims(table_name = table_name_chr_lp_mei, burn_in = 500)
cab_lp_mei_df <-
  run_qry_get_sims(table_name = table_name_cab_lp_mei, burn_in = 500)


# CABEZON
# use species specific vB params to create a weight at age data.frame that can be used to convert abundance at age to biomass at age
species_name_cab <- cab_nd_mei_nice[[2]]
cab_waa_g <- sim_species_parms[[species_name_cab]][["biom_const"]] * sim_species_derived_vars[[species_name_cab]][["length_at_age"]] ^
                    sim_species_parms[[species_name_cab]][["biom_exp"]]

cab_waa_g_df <- data.frame(age = 1: sim_species_parms[[species_name_cab]]$A_max,
                           waa_g = cab_waa_g)

# use functions to create useful data structures
# FLEP to F reference tables by species
spec_fleps_2_fs_cab <-
  get_fs_from_fleps(species = species_name,
                    flep_inds = species_flep_fished_ind,
                    f_vals =  fishing_mortality_values) 

# CHINA ROCKFISH
# use species specific vB params to create a weight at age data.frame that can be used to convert abundance at age to biomass at age
species_name_chr <- chr_nd_mei_nice[[2]]
chr_waa_g <- sim_species_parms[[species_name_chr]][["biom_const"]] * sim_species_derived_vars[[species_name_chr]][["length_at_age"]] ^
                    sim_species_parms[[species_name_chr]][["biom_exp"]]

chr_waa_g_df <- data.frame(age = 1: sim_species_parms[[species_name_chr]]$A_max,
                           waa_g = chr_waa_g)

# use functions to create useful data structures
# FLEP to F reference tables by species
spec_fleps_2_fs_chr <-
  get_fs_from_fleps(species = species_name_chr,
                    flep_inds = species_flep_fished_ind,
                    f_vals =  fishing_mortality_values) 

F_chr_nd_white_df_maxes_yield <- get_f_maxes_yield(input_df = chr_nd_white_df,
                                         fleps_2_fs = spec_fleps_2_fs_chr)

F_cab_nd_white_df_maxes_yield <- get_f_maxes_yield(input_df = cab_nd_white_df,
                                         fleps_2_fs = spec_fleps_2_fs_cab)

F_chr_lp_white_df_maxes_yield <- get_f_maxes_yield(input_df = chr_lp_white_df,
                                         fleps_2_fs = spec_fleps_2_fs_chr)

F_cab_lp_white_df_maxes_yield <- get_f_maxes_yield(input_df = cab_lp_white_df,
                                         fleps_2_fs = spec_fleps_2_fs_cab)

F_chr_nd_mei_df_maxes_yield <- get_f_maxes_yield(input_df = chr_nd_mei_df, 
                                                 fleps_2_fs = spec_fleps_2_fs_chr)

F_cab_nd_mei_df_maxes_yield <- get_f_maxes_yield(input_df = cab_nd_mei_df, 
                                                 fleps_2_fs = spec_fleps_2_fs_cab)

F_chr_lp_mei_df_maxes_yield <- get_f_maxes_yield(input_df = chr_lp_mei_df, 
                                                 fleps_2_fs = spec_fleps_2_fs_chr)

F_cab_lp_mei_df_maxes_yield <- get_f_maxes_yield(input_df = cab_lp_mei_df, 
                                                 fleps_2_fs = spec_fleps_2_fs_cab)


##----##
# No dispersal
# white
chr_nd_white_out_fs <- df_add_fs_out (input_df = chr_nd_white_df,
                         spec_fleps_2_fs = spec_fleps_2_fs_chr,
                         F_maxes_yield = F_chr_nd_white_df_maxes_yield, 
                         waa_g_df = chr_waa_g_df)

cab_nd_white_out_fs <- df_add_fs_out (input_df = cab_nd_white_df,
                         spec_fleps_2_fs = spec_fleps_2_fs_cab,
                         F_maxes_yield = F_cab_nd_white_df_maxes_yield, 
                         waa_g_df = cab_waa_g_df)
# enso
chr_nd_mei_out_fs <- df_add_fs_out (input_df = chr_nd_mei_df ,
                         spec_fleps_2_fs = spec_fleps_2_fs_chr,
                         F_maxes_yield = F_chr_nd_mei_df_maxes_yield, 
                         waa_g_df = chr_waa_g_df)

cab_nd_mei_out_fs <- df_add_fs_out (input_df = cab_nd_mei_df,
                         spec_fleps_2_fs = spec_fleps_2_fs_cab,
                         F_maxes_yield = F_cab_nd_mei_df_maxes_yield, 
                         waa_g_df = cab_waa_g_df)

# Larval pool
# white
chr_lp_white_out_fs <- df_add_fs_out (input_df = chr_lp_white_df,
                         spec_fleps_2_fs = spec_fleps_2_fs_chr,
                         F_maxes_yield = F_chr_lp_white_df_maxes_yield, 
                         waa_g_df = chr_waa_g_df)

cab_lp_white_out_fs <- df_add_fs_out (input_df = cab_lp_white_df,
                         spec_fleps_2_fs = spec_fleps_2_fs_cab,
                         F_maxes_yield = F_cab_lp_white_df_maxes_yield, 
                         waa_g_df = cab_waa_g_df)

# enso
chr_lp_mei_out_fs <- df_add_fs_out (input_df = chr_lp_mei_df ,
                         spec_fleps_2_fs = spec_fleps_2_fs_chr,
                         F_maxes_yield = F_chr_lp_mei_df_maxes_yield, 
                         waa_g_df = chr_waa_g_df)

cab_lp_mei_out_fs <- df_add_fs_out (input_df = cab_lp_mei_df,
                         spec_fleps_2_fs = spec_fleps_2_fs_cab,
                         F_maxes_yield = F_cab_lp_mei_df_maxes_yield, 
                         waa_g_df = cab_waa_g_df)

# White
# No dispersal
chr_nd_white_med_yield <- get_median_yield(input_df = chr_nd_white_out_fs) 

cab_nd_white_med_yield <- get_median_yield(input_df = cab_nd_white_out_fs) 

# larval pool dispersal 
chr_lp_white_med_yield <- get_median_yield(input_df = chr_lp_white_out_fs) 

cab_lp_white_med_yield <- get_median_yield(input_df = cab_lp_white_out_fs) 

# ENSO
# No dispersal
chr_nd_mei_med_yield <- get_median_yield(input_df = chr_nd_mei_out_fs) 

cab_nd_mei_med_yield <- get_median_yield(input_df = cab_nd_mei_out_fs) 

# larval pool dispersal
chr_lp_mei_med_yield <- get_median_yield(input_df = chr_lp_mei_out_fs) 

cab_lp_mei_med_yield <- get_median_yield(input_df = cab_lp_mei_out_fs) 

# White
# No dispersal
chr_nd_white_med_biomass <- get_median_biomass(input_df = chr_nd_white_out_fs) 

cab_nd_white_med_biomass <- get_median_biomass(input_df = cab_nd_white_out_fs) 

chr_lp_white_med_biomass <- get_median_biomass(input_df = chr_lp_white_out_fs) 

cab_lp_white_med_biomass <- get_median_biomass(input_df = cab_lp_white_out_fs) 

# ENSO
# No dispersal
chr_nd_mei_med_biomass <- get_median_biomass(input_df = chr_nd_mei_out_fs) 

cab_nd_mei_med_biomass <- get_median_biomass(input_df = cab_nd_mei_out_fs) 

# larval pool dispersal
chr_lp_mei_med_biomass <- get_median_biomass(input_df = chr_lp_mei_out_fs) 

cab_lp_mei_med_biomass <- get_median_biomass(input_df = cab_lp_mei_out_fs) 


chr_nd_white_out_yld_bm_sims_f_frac_yr <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = chr_nd_white_out_fs)
cab_nd_white_out_yld_bm_sims_f_frac_yr <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = cab_nd_white_out_fs)
chr_lp_white_out_yld_bm_sims_f_frac_yr <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = chr_lp_white_out_fs)
cab_lp_white_out_yld_bm_sims_f_frac_yr <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = cab_lp_white_out_fs)
chr_nd_mei_out_yld_bm_sims_f_frac_yr <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = chr_nd_mei_out_fs)
cab_nd_mei_out_yld_bm_sims_f_frac_yr <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = cab_nd_mei_out_fs)
chr_lp_mei_out_yld_bm_sims_f_frac_yr <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = chr_lp_mei_out_fs)
cab_lp_mei_out_yld_bm_sims_f_frac_yr <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = cab_lp_mei_out_fs)

##----##

chr_nd_white_sum_stats_f_frac_df <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = chr_nd_white_out_yld_bm_sims_f_frac_yr,
    median_yield = chr_nd_white_med_yield,
    median_biomass = chr_nd_white_med_biomass
  )

cab_nd_white_sum_stats_f_frac_df <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = cab_nd_white_out_yld_bm_sims_f_frac_yr,
    median_yield = cab_nd_white_med_yield,
    median_biomass = cab_nd_white_med_biomass
  )

chr_lp_white_sum_stats_f_frac_df <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = chr_lp_white_out_yld_bm_sims_f_frac_yr,
    median_yield = chr_lp_white_med_yield,
    median_biomass = chr_lp_white_med_biomass
  )

cab_lp_white_sum_stats_f_frac_df <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = cab_lp_white_out_yld_bm_sims_f_frac_yr,
    median_yield = cab_lp_white_med_yield,
    median_biomass = cab_lp_white_med_biomass
  )

chr_nd_mei_sum_stats_f_frac_df <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = chr_nd_mei_out_yld_bm_sims_f_frac_yr,
    median_yield = chr_nd_mei_med_yield,
    median_biomass =   chr_nd_mei_med_biomass
  )

cab_nd_mei_sum_stats_f_frac_df <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = cab_nd_mei_out_yld_bm_sims_f_frac_yr,
    median_yield = cab_nd_mei_med_yield,
    median_biomass =   cab_nd_mei_med_biomass
  )

chr_lp_mei_sum_stats_f_frac_df <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = chr_lp_mei_out_yld_bm_sims_f_frac_yr,
    median_yield = chr_lp_mei_med_yield,
    median_biomass =   chr_lp_mei_med_biomass
  )

cab_lp_mei_sum_stats_f_frac_df <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = cab_lp_mei_out_yld_bm_sims_f_frac_yr,
    median_yield = cab_lp_mei_med_yield,
    median_biomass =   cab_lp_mei_med_biomass
  )

Supp Fig 1

Biomass and yield contours for Blue rockfish

  • white noise

  • no dispersal and larval pool

# Species: Blue rockfish
# Dispersal: No Dispersal/Larval Pool
# Environmental noise: White

sim_res_out_brf_nd_white <-
  qry_get_sim_data_for_analysis(table_name = table_name_brf_nd_white)
sim_res_out_brf_lp_white <-
  qry_get_sim_data_for_analysis(table_name = table_name_brf_lp_white)


F_brf_maxes_yield_brf_nd_white <-
  get_f_maxes_yield(input_df = sim_res_out_brf_nd_white,
                    fleps_2_fs = spec_fleps_2_fs_brf)
F_brf_maxes_yield_brf_lp_white <-
  get_f_maxes_yield(input_df = sim_res_out_brf_lp_white,
                    fleps_2_fs = spec_fleps_2_fs_brf)

out_fs_brf_nd_white <- df_add_fs_out(
  input_df = sim_res_out_brf_nd_white,
  spec_fleps_2_fs = spec_fleps_2_fs_brf,
  F_maxes_yield = F_brf_maxes_yield_brf_nd_white,
  waa_g_df = brf_waa_g_df
)
out_fs_brf_lp_white <- df_add_fs_out(
  input_df = sim_res_out_brf_lp_white,
  spec_fleps_2_fs = spec_fleps_2_fs_brf,
  F_maxes_yield = F_brf_maxes_yield_brf_lp_white,
  waa_g_df = brf_waa_g_df
)

med_yield_brf <- get_median_yield(input_df = out_fs_brf_nd_white) 

med_biomass_brf <- get_median_biomass(input_df = out_fs_brf_nd_white)


out_yld_bm_sims_f_frac_yr_brf_nd_white <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = out_fs_brf_nd_white)
out_yld_bm_sims_f_frac_yr_brf_lp_white <-
  df_yld_bm_stats_by_sim_f_frac_year(input_df = out_fs_brf_lp_white)

sum_stats_f_frac_df_brf_nd_white <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = out_yld_bm_sims_f_frac_yr_brf_nd_white,
    median_yield = med_yield_brf,
    median_biomass = med_biomass_brf
  )
sum_stats_f_frac_df_brf_lp_white <-
  summary_stats_by_f_frac(
    out_sim_f_frac_yr = out_yld_bm_sims_f_frac_yr_brf_lp_white,
    median_yield = med_yield_brf,
    median_biomass = med_biomass_brf
  )

## Plotting
p1 <-
  plot_contour_bm_below_thresh(sum_stats_f_frac_df = sum_stats_f_frac_df_brf_nd_white,
                               spec_disp_noise = brf_nd_white_nice)

p2 <-
  plot_contour_yld_below_thresh(sum_stats_f_frac_df = sum_stats_f_frac_df_brf_nd_white,
                                spec_disp_noise = brf_nd_white_nice)
p3 <-
  plot_contour_bm_below_thresh(sum_stats_f_frac_df = sum_stats_f_frac_df_brf_lp_white,
                               spec_disp_noise = brf_lp_white_nice)

p4 <-
  plot_contour_yld_below_thresh(sum_stats_f_frac_df = sum_stats_f_frac_df_brf_lp_white,
                                spec_disp_noise = brf_lp_white_nice)
p_s1 <-
  plot_grid(
    p1,
    p2,
    p3,
    p4,
    labels = letters[1:4],
    label_size = 9,
    byrow = FALSE
  )

print(p_s1)

ggsave(filename = "output/supp_figs/sup_fig1-brf-biomass-yield-white-nd-lp.png", 
       plot = p_s1, 
       width = 12,
       height = 8)
ggsave(filename = "output/supp_figs/sup_fig1-brf-biomass-yield-white-nd-lp.svg", 
       plot = p_s1, 
       width = 12,
       height = 8)

# DBI::dbDisconnect(conn = con, shutdown=TRUE)

Supp Fig 2

China rockfish and Cabezon – Biomass and Yield contour plots

  • ENSO noise

  • No dispersal

## Contour plot for each species
### China rockfish
# Biomass

p1 <- chr_nd_mei_sum_stats_f_frac_df[[2]] %>%
  filter(as.numeric(reserve_frac) <= max_reserve_frac, f_fmsy <= max_f_fmsy) %>%
  ggplot(.,
       aes(x = as.numeric(reserve_frac), y = f_fmsy, z = yrs_below_biomass_thresh)) +
  metR::geom_contour_fill(bins = 12) +
  ggtitle(paste0(
        chr_nd_mei_nice[[2]],
        " ",
        chr_nd_mei_nice[[1]],
        " ",
        chr_nd_mei_nice[[3]],
        ":\nTime below mean median biomass "
      )) +
      xlab("Fraction of coastline in reserves") +
      ylab("Harvest Rate (F/Fmsy)") +
      scale_fill_viridis_c(option = "D", direction = -1) +
      geom_segment(aes(
        x = 0,
        xend = max_reserve_frac,
        y = 1,
        yend = 1
      ), colour = "black",
      linetype = 2) +
      scale_x_continuous(expand = expansion(mult = 0, add = 0), limits = c(0, 0.3)) +
      scale_y_continuous(expand = expansion(mult = 0, add = 0)) +
      theme_bw()

# Yield
p2 <- chr_nd_mei_sum_stats_f_frac_df[[2]] %>%
  filter(as.numeric(reserve_frac) <= max_reserve_frac, f_fmsy <= max_f_fmsy) %>%
  ggplot(.,
       aes(x = as.numeric(reserve_frac), y = f_fmsy, z = yrs_below_yield_thresh)) +
  metR::geom_contour_fill(bins = 12) + 
  ggtitle(paste0(
        chr_nd_mei_nice[[2]],
        " ",
        chr_nd_mei_nice[[1]],
        " ",
        chr_nd_mei_nice[[3]],
        ":\nTime below median max yield"
      )) +
      xlab("Fraction of coastline in reserves") +
      ylab("Harvest Rate (F/Fmsy)") +
      scale_fill_viridis_c(option = "D", direction = -1) +
      geom_segment(aes(
        x = 0,
        xend = max_reserve_frac,
        y = 1,
        yend = 1
      ), colour = "black",
      linetype = 2) +
      scale_x_continuous(expand = expansion(mult = 0, add = 0), limits = c(0, 0.3)) +
      scale_y_continuous(expand = expansion(mult = 0, add = 0)) +
      theme_bw()

### Cabezon
# Biomass
p3 <- cab_nd_mei_sum_stats_f_frac_df[[2]] %>%
  filter(as.numeric(reserve_frac) <= max_reserve_frac, f_fmsy <= max_f_fmsy) %>%
  ggplot(.,
       aes(x = as.numeric(reserve_frac), y = f_fmsy, z = yrs_below_biomass_thresh)) +
  metR::geom_contour_fill(bins = 12) +
  ggtitle(paste0(
        cab_nd_mei_nice[[2]],
        " ",
        cab_nd_mei_nice[[1]],
        " ",
        cab_nd_mei_nice[[3]],
        ":\nTime below mean median biomass "
      )) +
      xlab("Fraction of coastline in reserves") +
      ylab("Harvest Rate (F/Fmsy)") +
      scale_fill_viridis_c(option = "D", direction = -1) +
      geom_segment(aes(
        x = 0,
        xend = max_reserve_frac,
        y = 1,
        yend = 1
      ), colour = "black",
      linetype = 2) +
      scale_x_continuous(expand = expansion(mult = 0, add = 0), limits = c(0, 0.3)) +
      scale_y_continuous(expand = expansion(mult = 0, add = 0)) +
      theme_bw()

# Yield
p4 <- cab_nd_mei_sum_stats_f_frac_df[[2]] %>%
  filter(as.numeric(reserve_frac) <= max_reserve_frac, f_fmsy <= max_f_fmsy) %>%
  ggplot(.,
       aes(x = as.numeric(reserve_frac), y = f_fmsy, z = yrs_below_yield_thresh)) +
  metR::geom_contour_fill(bins = 12) + 
  ggtitle(paste0(
        cab_nd_mei_nice[[2]],
        " ",
        cab_nd_mei_nice[[1]],
        " ",
        cab_nd_mei_nice[[3]],
        ":\nTime below median max yield"
      )) +
      xlab("Fraction of coastline in reserves") +
      ylab("Harvest Rate (F/Fmsy)") +
      scale_fill_viridis_c(option = "D", direction = -1) +
      geom_segment(aes(
        x = 0,
        xend = max_reserve_frac,
        y = 1,
        yend = 1
      ), colour = "black",
      linetype = 2) +
      scale_x_continuous(expand = expansion(mult = 0, add = 0), limits = c(0, 0.3)) +
      scale_y_continuous(expand = expansion(mult = 0, add = 0)) +
      theme_bw()

p_s2 <- plot_grid(p1, p2, p3, p4, labels = letters[1:4], label_size = 12, byrow = FALSE)

print(p_s2)

ggsave(filename = "output/supp_figs/sup-fig2-cab-chr-nodisp-mei-yld-biom.png", 
       plot = p_s2, 
       width = 12,
       height = 8)
ggsave(filename = "output/supp_figs/sup-fig2-cab-chr-nodisp-mei-yld-biom.svg", 
       plot = p_s2, 
       width = 12,
       height = 8)

Supp Fig 3

China rockfish and Cabezon – Biomass and Yield contour plots

  • ENSO noise

  • larval pool dispersal

## Contour plot for each species
### China rockfish
# Biomass

p1 <- chr_lp_mei_sum_stats_f_frac_df[[2]] %>%
  filter(as.numeric(reserve_frac) <= max_reserve_frac, f_fmsy <= max_f_fmsy) %>%
  ggplot(.,
       aes(x = as.numeric(reserve_frac), y = f_fmsy, z = yrs_below_biomass_thresh)) +
  metR::geom_contour_fill(bins = 12) +
  ggtitle(paste0(
        chr_lp_mei_nice[[2]],
        " ",
        chr_lp_mei_nice[[1]],
        " ",
        chr_lp_mei_nice[[3]],
        ":\nTime below mean median biomass "
      )) +
      xlab("Fraction of coastline in reserves") +
      ylab("Harvest Rate (F/Fmsy)") +
      scale_fill_viridis_c(option = "D", direction = -1) +
      geom_segment(aes(
        x = 0,
        xend = max_reserve_frac,
        y = 1,
        yend = 1
      ), colour = "black",
      linetype = 2) +
      scale_x_continuous(expand = expansion(mult = 0, add = 0), limits = c(0, 0.3)) +
      scale_y_continuous(expand = expansion(mult = 0, add = 0)) +
      theme_bw()

# Yield
p2 <- chr_lp_mei_sum_stats_f_frac_df[[2]] %>%
  filter(as.numeric(reserve_frac) <= max_reserve_frac, f_fmsy <= max_f_fmsy) %>%
  ggplot(.,
       aes(x = as.numeric(reserve_frac), y = f_fmsy, z = yrs_below_yield_thresh)) +
  metR::geom_contour_fill(bins = 12) + 
  ggtitle(paste0(
        chr_lp_mei_nice[[2]],
        " ",
        chr_lp_mei_nice[[1]],
        " ",
        chr_lp_mei_nice[[3]],
        ":\nTime below median max yield"
      )) +
      xlab("Fraction of coastline in reserves") +
      ylab("Harvest Rate (F/Fmsy)") +
      scale_fill_viridis_c(option = "D", direction = -1) +
      geom_segment(aes(
        x = 0,
        xend = max_reserve_frac,
        y = 1,
        yend = 1
      ), colour = "black",
      linetype = 2) +
      scale_x_continuous(expand = expansion(mult = 0, add = 0), limits = c(0, 0.3)) +
      scale_y_continuous(expand = expansion(mult = 0, add = 0)) +
      theme_bw()

print(p1)

print(p2)

### Cabezon
# Biomass
p3 <- cab_lp_mei_sum_stats_f_frac_df[[2]] %>%
  filter(as.numeric(reserve_frac) <= max_reserve_frac, f_fmsy <= max_f_fmsy) %>%
  ggplot(.,
       aes(x = as.numeric(reserve_frac), y = f_fmsy, z = yrs_below_biomass_thresh)) +
  metR::geom_contour_fill(bins = 12) +
  ggtitle(paste0(
        cab_lp_mei_nice[[2]],
        " ",
        cab_lp_mei_nice[[1]],
        " ",
        cab_lp_mei_nice[[3]],
        ":\nTime below mean median biomass "
      )) +
      xlab("Fraction of coastline in reserves") +
      ylab("Harvest Rate (F/Fmsy)") +
      scale_fill_viridis_c(option = "D", direction = -1) +
      geom_segment(aes(
        x = 0,
        xend = max_reserve_frac,
        y = 1,
        yend = 1
      ), colour = "black",
      linetype = 2) +
      scale_x_continuous(expand = expansion(mult = 0, add = 0), limits = c(0, 0.3)) +
      scale_y_continuous(expand = expansion(mult = 0, add = 0)) +
      theme_bw()

# Yield
p4 <- cab_lp_mei_sum_stats_f_frac_df[[2]] %>%
  filter(as.numeric(reserve_frac) <= max_reserve_frac, f_fmsy <= max_f_fmsy) %>%
  ggplot(.,
       aes(x = as.numeric(reserve_frac), y = f_fmsy, z = yrs_below_yield_thresh)) +
  metR::geom_contour_fill(bins = 12) + 
  ggtitle(paste0(
        cab_lp_mei_nice[[2]],
        " ",
        cab_lp_mei_nice[[1]],
        " ",
        cab_lp_mei_nice[[3]],
        ":\nTime below median max yield"
      )) +
      xlab("Fraction of coastline in reserves") +
      ylab("Harvest Rate (F/Fmsy)") +
      scale_fill_viridis_c(option = "D", direction = -1) +
      geom_segment(aes(
        x = 0,
        xend = max_reserve_frac,
        y = 1,
        yend = 1
      ), colour = "black",
      linetype = 2) +
      scale_x_continuous(expand = expansion(mult = 0, add = 0), limits = c(0, 0.3)) +
      scale_y_continuous(expand = expansion(mult = 0, add = 0)) +
      theme_bw()

p_s3 <- plot_grid(p1, p2, p3, p4, labels = letters[1:4], label_size = 12, byrow = FALSE)

print(p_s3)

ggsave(filename = "output/supp_figs/sup-fig3-cab-chr-lp-disp-mei-yld-biom.png", 
       plot = p_s3, 
       width = 12,
       height = 8)
ggsave(filename = "output/supp_figs/sup-fig3-cab-chr-lp-disp-mei-yld-biom.svg", 
       plot = p_s3, 
       width = 12,
       height = 8)